The problem of finding the intersection of two circles in a plane is routine. At least, its easy if we are using numbers and even easier if we are using numbers and a computer algebra system (CAS). However, if we just set up the symbolic problem and ask the CAS to give a general symbolic solution, it is horrendous. For example, let the first circle be centered at $(h,k)$ and have radius $r$. Then the locus of points is described by the equation $(x-h)^{2}+(y-k)^{2}=r^{2}$. Now we can have a second circle centered at $(m,n)$ with radius $s$. So our two equations are $$(x-h)^{2}+(y-k)^{2}=r^{2}$$ $$(x-m)^{2}+(y-n)^{2}=s^{2}$$ If the circles overlap such that they have two intersection points, then we can solve them simultaneously to get values for $(x_{1},y_{1})$ and $(x_{2},y_{2})$. Geogebra's CAS would not do this problem unless most of the constants were given numerical values. SAGE would do it, but the answer is not presentable in polite society. Wolfram Alpha, the online Mathematica tool, didn't choke, but also didn't give a useful result. So, what we will do in this section is show a viable way to approach this problem other than just getting a numerical computer solution. I am going to show three ways to solve this problem: a) The easy intersection, b) Intersect by Rotation Translation, c) Intersect using vectors. Method āaā, is absolutely essential to understanding the other two. Method ācā requires understanding a bit about vector projection, but is unquestionably the best way to go.
We will begin by putting one of the two circles at the origin and putting the other circle on the $x$-axis. Now, circle $A$ in figure 1 has the equation $x^{2}+y^{2}=r_{1}^{2}$ and circle $B$ has the equation $(x-d)^{2}+y^{2}=r_{2}^{2}$ where $d$ is the distance between centers. We will use the first equation to find $y^{2}=r_{1}^{2}-x^{2}$ and substitute that into the second equation. $$(x-d)^{2}+r_{1}^{2}-x^{2}=r_{2}^{2} \tag{1} \label{1}$$ Now expand $\eqref{1}$. $$x^{2}-2xd+d^{2}+r_{1}^{2}-x^{2}=r_{2}^{2}$$ Simplify to lose the $x^2$ term. $$-2xd+d^{2}+r_{1}^{2}=r_{2}^{2}$$ Solve for $x$. Remember, $x$ is $P_{x}$ which is one of the coordinate values of the intersection! $$x=\frac{d^{2}+r_{1}^{2}-r_{2}^{2}}{2d}$$ To get $y$, we can solve the equation of circle $A$ for $y$. $$y=\pm\sqrt{r_{1}^{2}-x^{2}}$$ Thus the intersections are $P,Q=(x,\pm y)$.
The general circles are not at the origin and on the axis. From here, we are going to define the general circle intersection points as operations rather than equations. The idea is that we will start with two circles on the plane, called $A$ and $B$. We will translate these so that circle $A$ is at the origin and then rotate them such that circle $B$ is on the axis. Then we will solve for the intersection just as in The Easy Intersection and then we will un-rotate and un-translate in order to put the intersections back where they belong.
Since these translations and rotations have nothing whatever to do with the inputs, $d,\,r_{1}$ and $r_{2}$, we can begin by defining these constants and plan to calculate $P$ and $Q$ just as in the prior section. $d=|B-A|$. $A$ is the center of circle $A$ and $B$ is the center of circle $B$.
Summary Plan
We have to understand how the translation and rotation are performed. Let's start with the translation, since it is easiest and has to be done first anyway. To both point, $A$ and $B$ we add $(-A_{x},-A_{y})$. $$Translation\;of\;A=\left(\begin{array}{c} A_{x}\\ A_{y} \end{array}\right)+\left(\begin{array}{c} -A_{x}\\ -A_{y} \end{array}\right)=\left(\begin{array}{c} 0\\ 0 \end{array}\right)$$ $$Translation\;of\;B=B'=\left(\begin{array}{c} B_{x}\\ B_{y} \end{array}\right)+\left(\begin{array}{c} -A_{x}\\ -A_{y} \end{array}\right)=\left(\begin{array}{c} B_{x}-A_{x}\\ B_{y}-A_{y} \end{array}\right)$$
Note that $B^{'}$ is not yet on the $x$-axis. It is translated, but not rotated. After the translation, in order to rotate point $B$ onto the axis, we need to know the angle $\theta$ that is between $\overline{AB}$ and the axis. To find angle $\theta$, since it is derived from a vector from the origin, we can use $\theta=atan2(B_{y},B_{x})$.
To rotate $B^{'}$ onto the $x$-axis multiply this clockwise rotation matrix time $B^{'}$. $$B^{''}=\left[\begin{array}{cc}\cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{array} \right] \left(\begin{array}{c} B_{x}^{'}\\ B_{y}^{'} \end{array}\right)$$ Next calculate intersection points $P$ and $Q$. This should be easy since one of the circle centers, $A^{'}$ is at the origin and the other one, $B^{''}$ is on the $x$-axis.
If $\theta$ is how much we rotate down, then $\theta$ is also how much we need to rotate back up again. To rotate the point back into position, we multiply by a counter-clockwise rotation matrix. (Remember to do both $P$ and $Q$). $$\left(\begin{array}{cc} cos\theta & -sin\theta\\ sin\theta & cos\theta \end{array}\right)\left(\begin{array}{c} P_{x}\\ P_{y} \end{array}\right)=P'$$ $$\left(\begin{array}{cc} cos\theta & -sin\theta\\ sin\theta & cos\theta \end{array}\right)\left(\begin{array}{c} Q_{x}\\ Q_{y} \end{array}\right)=Q'$$
Now we apply the un-translate transformation to get the original intersections.
Un-translate $P$.
$$\left(\begin{array}{c}
P_{x}^{'}\\
P_{y}^{'}
\end{array}\right)+\left(\begin{array}{c}
A_{x}\\
A_{y}
\end{array}\right)=\left(\begin{array}{c}
P_{x}^{'}+A_{x}\\
P_{y}^{'}+A_{y}
\end{array}\right)=\left(\begin{array}{c}
P_{x}^{''}\\
P_{y}^{''}
\end{array}\right)$$
Un-translate $Q$.
$$\left(\begin{array}{c}
Q_{x}^{'}\\
Q_{y}^{'}
\end{array}\right)+\left(\begin{array}{c}
A_{x}\\
A_{y}
\end{array}\right)=\left(\begin{array}{c}
Q_{x}^{'}+A_{x}\\
Q_{y}^{'}+A_{y}
\end{array}\right)=\left(\begin{array}{c}
Q_{x}^{''}\\
Q_{y}^{''}
\end{array}\right)$$
The example is a live Geogebra sheet to show that this works in all sectors. You can drag points $A$ and $B$ and the intersection will continue to be computed correctly.
You have read and understood the easy intersection . Therefore, we already know that our intersection points, $P$ and $Q$ have the following formulas. $$P=(x_{1},y_{1})=\left(\begin{array}{c} {\displaystyle \frac{d^{2}+r_{1}^{2}-r_{2}^{2}}{2d}}\\ \\ \sqrt{r_{1}^{2}-x_{1}^{2}} \end{array}\right)\qquad Q=(x_{2},y_{2})=\left(\begin{array}{c} {\displaystyle \frac{d^{2}+r_{1}^{2}-r_{2}^{2}}{2d}}\\ \\ -\sqrt{r_{1}^{2}-x_{1}^{2}} \end{array}\right)$$ We also know that $d=|B-A|$ which is the distance from point $A$ to point $B$. If we were to construct a line between points $P$ and $Q$, then the midpoint of that line, $M$, would be $$M=\left(\frac{d^{2}+r_{1}^{2}-r_{2}^{2}}{2d},\;0\right)$$ If $A$ was not at the origin and $B$ was not on the $x$-axis, then $|M|$ would be the distance along the vector from $A$ to $B$. Pretty clearly, the absolute value of $M$ is just the calculation $x_{1}$. In order to project along $A\rightarrow B$, we need the unit vector between $A$ and $B$. Call it $\mathbf{u}$. (Note: we have already computed the denominator of $\mathbf{u}$ and it is distance, $d$.) $$\mathbf{u}=\frac{(B-A)}{\left|(B-A)\right|}=\frac{\left(\begin{array}{c} B_{x}\\ B_{y} \end{array}\right)-\left(\begin{array}{c} A_{x}\\ A_{y} \end{array}\right)}{\sqrt{\left(B_{x}-A_{x}\right)^{2}-\left(B_{y}-A_{y}\right)^{2}}}$$ So when not along an axis, the coordinates of $M$ are given by: $$\left(\begin{array}{c} M_{x}\\ M_{y} \end{array}\right)=\left(\begin{array}{c} A_{x}\\ A_{y} \end{array}\right)+x_{1}\left(\begin{array}{c} u_{x}\\ u_{y} \end{array}\right).$$ But we are actually looking for $P$ and $Q$. The vector that is perpendicular to $\mathbf{u}$, call it $\mathbf{v}$, is given by $$\mathbf{v}=\left(\begin{array}{c} -\mathbf{u_{y}}\\ \mathbf{u_{x}} \end{array}\right)$$ Now we have the direction from point $M$ to point $P$. But, we also know that the distance in that direction if given by $y_{1}$ or if we wanted to go to $Q$, then $y_{2}$. We know both of these values from the Easy Intersection Thus: $$P=\left(\begin{array}{c} M_{x}\\ M_{y} \end{array}\right)+y_{1}\left(\begin{array}{c} \mathbf{v_{x}}\\ \mathbf{v_{y}} \end{array}\right)\qquad Q=\left(\begin{array}{c} M_{x}\\ M_{y} \end{array}\right)-y_{1}\left(\begin{array}{c} \mathbf{v_{x}}\\ \mathbf{v_{y}} \end{array}\right)$$